/*
 * CorrCoeff.h
 *
 *  Created on: Jan 9, 2012
 *      Author: dblyth
 */

#ifndef CORRCOEFF_H_
#define CORRCOEFF_H_

template<class T>
class CorrCoeff {
public:
	CorrCoeff(double* ref, double threshold);
	virtual ~CorrCoeff();

	void
	getCoeffsOverThres(T* stream, int streamSize, bool* isOver, double* coeff = 0);
	bool hasCoeffOverThres(T* stream, int streamSize);

private:
	double* ref;
	int refSize;
	double refSum;
	double ref2Sum;
	double beta;
	double gamma;

};

#include <cmath>

template<class T>
CorrCoeff<T>::CorrCoeff(double* ref, double threshold) {
	refSum = 0;
	refSize = 21;

	for (int i = 0; i < refSize; i++)
		refSum += ref[i];

	double factor = 1.0 / refSum;
	this->ref = new double[refSize];

	for (int i = 0; i < refSize; i++) {
		this->ref[i] = factor * ref[i];
		ref2Sum += this->ref[i] * this->ref[i];
	}

	refSum = 1.0;

	gamma = refSize * ref2Sum - 1;
	beta = gamma * threshold * threshold;

	for (int i = 0; i < refSize; i++)
		this->ref[i] *= refSize;
}

template<class T>
CorrCoeff<T>::~CorrCoeff() {
	delete[] ref;
}

template<class T>
void CorrCoeff<T>::getCoeffsOverThres(T* stream, int streamSize, bool* isOver, double* coeff) {
	if (streamSize < refSize)
		return;

	double dataSum = 0;
	double data2Sum = 0;
	double dataRefSum;

	for (int i = 0; i < refSize; i++) {
		dataSum += stream[i];
		data2Sum += stream[i] * stream[i];
	}

	unsigned int const NSMRS = streamSize - refSize;

	for (int i = 0; i <= NSMRS; i++) {
		T* const tempStream = &stream[i];

		double DRSTRS = tempStream[0] * ref[0] + tempStream[1] * ref[1] + tempStream[2] * ref[2] + tempStream[3]
				* ref[3] + tempStream[4] * ref[4] + tempStream[5] * ref[5] + tempStream[6] * ref[6] + tempStream[7]
				* ref[7] + tempStream[8] * ref[8] + tempStream[9] * ref[9] + tempStream[10] * ref[10] + tempStream[11]
				* ref[11] + tempStream[12] * ref[12] + tempStream[13] * ref[13] + tempStream[14] * ref[14]
				+ tempStream[15] * ref[15] + tempStream[16] * ref[16] + tempStream[17] * ref[17] + tempStream[18]
				* ref[18] + tempStream[19] * ref[19] + tempStream[20] * ref[20];

		if (DRSTRS > dataSum) {
			double DRSTRSMDS = DRSTRS - dataSum;
			double RSTD2SMDS2 = refSize * data2Sum - dataSum * dataSum;
			if (DRSTRSMDS * DRSTRSMDS > beta * RSTD2SMDS2) {
				if (isOver)
					isOver[i] = true;
				if (coeff) {
					coeff[i] = DRSTRSMDS / pow(gamma * RSTD2SMDS2, 0.5);
				}
			} else {
				if (isOver)
					isOver[i] = false;
				if (coeff)
					coeff[i] = 0.;
			}
		} else {
			if (isOver)
				isOver[i] = false;
			if (coeff)
				coeff[i] = 0.;
		}

		if (i < NSMRS) {
			dataSum += tempStream[refSize] - tempStream[0];
			data2Sum += tempStream[refSize] * tempStream[refSize] - tempStream[0] * tempStream[0];
		}
	}
}

template<class T>
bool CorrCoeff<T>::hasCoeffOverThres(T* stream, int streamSize) {
	if (streamSize < refSize)
		return false;

	double dataSum = 0;
	double data2Sum = 0;
	double dataRefSum;

	for (int i = 0; i < refSize; i++) {
		dataSum += stream[i];
		data2Sum += stream[i] * stream[i];
	}

	unsigned int const NSMRS = streamSize - refSize;

	for (int i = 0; i <= NSMRS; i++) {
		T* const tempStream = &stream[i];

		double DRSTRS = tempStream[0] * ref[0] + tempStream[1] * ref[1] + tempStream[2] * ref[2] + tempStream[3]
				* ref[3] + tempStream[4] * ref[4] + tempStream[5] * ref[5] + tempStream[6] * ref[6] + tempStream[7]
				* ref[7] + tempStream[8] * ref[8] + tempStream[9] * ref[9] + tempStream[10] * ref[10] + tempStream[11]
				* ref[11] + tempStream[12] * ref[12] + tempStream[13] * ref[13] + tempStream[14] * ref[14]
				+ tempStream[15] * ref[15] + tempStream[16] * ref[16] + tempStream[17] * ref[17] + tempStream[18]
				* ref[18] + tempStream[19] * ref[19] + tempStream[20] * ref[20];

		if (DRSTRS > dataSum) {
			double DRSTRSMDS = DRSTRS - dataSum;
			double RSTD2SMDS2 = refSize * data2Sum - dataSum * dataSum;
			if (DRSTRSMDS * DRSTRSMDS > beta * RSTD2SMDS2) {
				return true;
			}
		}

		if (i < NSMRS) {
			dataSum += tempStream[refSize] - tempStream[0];
			data2Sum += tempStream[refSize] * tempStream[refSize] - tempStream[0] * tempStream[0];
		}
	}

	return false;
}

#endif /* CORRCOEFF_H_ */
